function [U, S, V, out] = MMBS( D, lambda, para )

model.lambda = 2*lambda;
params.tol = para.tol;

%% Setup
d1 = size(D, 1);
d2 = size(D, 2);

[row, col, D] = find(D);

% Options
compute_predictions = true;

%% Optimization algorithm is the trust-region algorithm
optim_algo = @tr_matrix_completion;

%% Paramters
params.vtol = 1e-6;
params.grad_tol = 1e-4;
params.dg_tol = 1e-3;
params.dg_vtol = 1e-4;
params.smax_tol = 1e-3;
params.max_iter = 10000;
params.max_iter_tr = 50;
params.verb = true; % Show output
params.compute_predictions = compute_predictions; % Compute predictions each iteration if asked
params = default_params(params);
 
% data 
data.ls.rows = row;
data.ls.cols = col;
data.ls.entries = D;
data.ls.nentries = length(D);
data.ts.rows = para.test.col;
data.ts.cols = para.test.row;
data.ts.entries = para.test.data;
data.ts.nentries = length(para.test.data);

data.nentries = length(D);
data.d1 = d1;
data.d2 = d2;

%% Creating a sparse structure
sparse_structure  = sparse(row, col, length(row), d1, d2);
sparse_structure2 = sparse(row, col, length(row), d1, d2);

%% Setup initial model
model.d1 = d1;
model.d2 = d2;
model.p = 1; 
model.pmax = min(d1, d2);
model.sparse_structure = sparse_structure;
model.sparse_structure2 = sparse_structure2;

%% Random Initialization
G = randn(model.d1, model.p); H = randn(model.d2, model.p);
[Qg, Rg] = qr(G, 0);
[Qh, Rh] = qr(H, 0);
[q1, b, q2] = svd(Rg * Rh');

model.U = Qg * q1;
model.V = Qh * q2;
model.B = b;

%% Run algorithm
Time = tic;
[model, info] = low_rank_optimization(optim_algo, data, model, params);
Time = toc(Time);

U = model.U;
V = model.V;
S = model.B;

out.obj = info.costs;
out.RMSE = sqrt(info.test_error/length(para.test.data));
out.Time = info.iter_time;

Time = Time/length(out.obj);
Time = Time*(1:length(out.obj));
out.Time = Time;
out.rank = size(S, 1);

end
